Theoretical description of phase coexistence in model Ceo 
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We have investigated the phase diagram of a pair interaction model of Ceo fullerene [L. A. Girifalco, 
J. Phys. Chem. 96, 858 (1992)], in the framework provided by two integral equation theories of the 
liquid state, namely the Modified Hypernetted Chain (MHNC) implemented under a global thermo- 
dynamic consistency constraint, and the Self-Consistent Ornstein-Zernike Approximation (SCOZA), 
and by a Perturbation Theory (PT) with various degrees of refinement, for the free energy of the 
solid phase. 

We present an extended assessment of such theories as set against a recent Monte Carlo study of 
the same model [D. Costa, G. Pellicane, C. Caccamo, and M. C. Abramo, J. Chem. Phys. 118, 304 
(2003)]. We have compared the theoretical predictions with the corresponding simulation results 
for several thermodynamic properties like the free energy, the pressure, and the internal energy. 
Then we have determined the phase diagram of the model, by using either the SCOZA, or the 
MHNC, or the PT predictions for one of the coexisting phases, and the simulation data for the 
other phase, in order to separately ascertain the accuracy of each theory. It turns out that the 
overall appearance of the phase portrait is reproduced fairly well by all theories, with remarkable 
accuracy as for the melting line and the solid-vapor equilibrium. All theories show a more or less 
pronounced discrepancy with the simulated fluid-solid coexistence pressure, above the triple point. 
The MHNC and SCOZA results for the liquid-vapor coexistence, as well as for the corresponding 
critical points, are quite accurate; the SCOZA tends to underestimate the density corresponding to 
the freezing line. All results are discussed in terms of the basic assumptions underlying each theory. 

We have selected the MHNC for the fluid and the first-order PT for the solid phase, as the 
most accurate tools to investigate the phase behavior of the model in terms of purely theoretical 
approaches. It emerges that the use of different procedures to characterize the fluid and the solid 
phases provides a semiquantitative reproduction of the thermodynamic properties of the Ceo model 
at issue. 

The overall results appear as a robust benchmark for further theoretical investigations on higher 
order Cn>60 fuUerenes, as well as on other fuUerene-related materials, whose description can be 
based on a modelization similar to that adopted in this work. 

PACS numbers: 61.20.Gy, 61.48.+C, 64.70.-p 



I. INTRODUCTION 

A current model for Ceo is based on a "smeared out" 
spherical representation of the fullerene molecules, early 
proposed by Girifalco 0, where the carbon atoms give 
rise to a uniform interaction distributed over the molec- 
ular cage surface. A straightforward integration leads to 
an analytical central pair potential, characterized by a 
harsh repulsive core followed by a deep attractive well, 
which rapidly decays with the interparticle distance. 

A detailed characterization of the thermodynamic 
properties of the Girifalco model is of relevant interest 
in several respects: for instance, first-principle studies of 
the Cgo interaction give results very similar to those pre- 
dicted through the Girifalco model ^ |3|| ; the latter has 
been widely used to model C„>6o fuUerenes with n = 70, 
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76, and 84 (see 0,|EIEI3 and references cited therein), as 
well as other hollow nanoparticles like carbon onions and 
metal dichalcogenides like GaAs and CdSe (also termed 
inorganic fuUerenes) P| . Umiguchi and coworkers 9] pro- 
posed a semiempirical modification of the Girifalco inter- 
action to account for the behavior of solid Cgo at low tem- 
perature; such a generalization might prove useful for the 
analysis of the lattice properties under impurity doping. 

As far as a more speculative analysis of the Girifalco 
model is concerned, several intriguing aspects are related 
to the overall appearance of its phase portrait, character- 
ized by a narrow liquid pocket that extends only a few 
tens degrees over the triple point temperature |0, 0| . 
This borderline behavior, as for the existence of a stable 
liquid phase, basically depends on the interplay between 
entropic demands, imposed by excluded-volume effects 
at short distance, and energetic contributions, due to the 
attractive part of the interaction potential. We may re- 
call, as an indication of the implied subtleties, the early 
controversy on the location of the liquid-vapor critical 
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point 10, , clarified in successive studies 0, 0| , and 
the discrepancies on the position of the freezing Hne and 
of the triple point, related to the procedure ado pted for 
the calculation of the liquid-solid equilibrium . 
Recently (and unexpectedly, on the basis of a previ- 
ous study on the Lennard- Jones fluid |E3), Fartaria and 
coworkers have found that even the same simulation 
approach leads to distinct predictions on the freezing line 
of the Girifalco model, provided that different starting 
thermodynamic conditions are employed. 

We have recently reconciled such disparate results in 
the context of extensive Monte Carlo calculations of the 
free energies of both the fluid and the solid phases of the 
model [ill ITsI l . We have concluded that the solid-fluid 
equilibrium is strongly affected by the deep, short-range 
attractive well in the interaction potential, at variance 
with systems whose freezing behavior is essentially dom- 
inated by steric effects. As a consequence, the freezing 
transition of the fluid is driven to lower densities (and 
the liquid pocket is reduced accordingly), mainly by en- 
ergetic effects, in agreement with a common scenario re- 
cently proposed for fluids interacting through short-range 
forces [l9| . 

The unusual aspects of the phase behavior have oc- 
casioned several studies on the Cgo model by means of 
refined theoretical tools, like for instance integral equa- 
tion theories for the fluid phase flsL ITp . various density 
functional approximations [l^ l2Ct l2l|. and the hierar- 
chical reference theory p^ . Preliminary studies have 
been recently carried out, based on the Modified Hy- 
pernetted Chain approach (MHNC, 23]), solved jinder 
a global thermodynamic consistency constraint 24] , and 
on the Self-Consistent Ornstein-Zernike Approximation 
(SCOZA m m 113) m. However, such investiga- 
tions were often compared to simulation results later re- 
vised; for example, several authors 13, 16,, 2§] have esti- 
mated the freezing conditions on the basis of some well- 
known one-phase structural indicators, as the Hansen- 
Verlet 29], or the Giaquinta and Giunta js^] criteria, 
whose limited applicability in the context of "energetic" 
fluids has been discussed in QA]. 

Hinging on the simulation study of Ref. |llj ] , it is now 
worth reconsidering a detailed analysis of theoretical pre- 
dictions for the phase diagram of the envisaged model. 
Following our preliminary investigations [2^. |2^ we shall 
adopt the MHNC and SCOZA theories to characterize 
the thermodynamic properties of the fluid phase. As far 
as the free energy of the solid phase is concerned, we 
use a Perturbation Theory (PT), discussing several de- 
grees of refinement, as fully described in the text. Recent 
studies (see 0, Is^] and references) have demonstrated 
that the PT accurately describes the solid phase and the 
solid-fluid transition of models characterized by short- 
range interactions, as the depletion potentials resulting 
from the effective one-component representation of hard- 
sphere mixtures. A second-order expansion for the PT, 
also analyzed in this work, has been used in Ref. jS^j to 
investigate the phase behavior of several simple models 



for globular protein solutions. 

This paper is organized as follows: the theoretical ap- 
proaches are described in section II. Results are reported 
and discussed in Section HI. A short overview of our re- 
sults and the conclusions are drawn in Section IV. 



II. MODEL AND THEORETICAL 
APPROACHES 



A. The Girifalco model potential 

The Girifalco potential is well known in the fuUerene 
literature, so we recall only its analytical expression 
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where s = r/d, ai = N'^A/12d^, and a2 = N'^B/90d^'^; 
N and d are the number of carbon atoms and the di- 
ameter, respectively, of the fuUerene particles, A = 32 x 
10^^" ergcm^ and B = 55.77 x 10^^°^ ergcm^^ are con- 
stants entering the Lennard- Jones 12-6 potential through 
which two carbon sites on different spherical molecules 
are assumed to interact m. For Ceo, d = 0.71 nm while 
the node of the potential (Q) , the minimum, and its po- 
sition, are tq ~ 0.959 nm, e ~ 0.444 x 10"^^ erg and 
Tmin = 1.005 nm, respectively. 



B. Perturbation theory of the solid phase 

The perturbation approach is based on a separation 
of the potential Q into a reference (purely repulsive) 
part, tirof(''), and a residual, attractive part, t;pcrt(?')- 
The reference system is then approximated by a solid of 
hard spheres with the same crystallographic structure, by 
means of a suitable definition of the hard-core diameter. 
Given these positions, the free energy of the system can 
be expanded around the free energy of the hard-sphere 
crystal, F^s, to obtain at first order [s^ ]: 
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"port(r) ffhsWdr, (2) 



where TV is the number of particles, (3 the inverse of the 
temperature T in units of the Boltzmann constant, and p 
the number density; the second term is the thermal aver- 
age of the perturbation energy over the reference system, 
where gj^g (r) is the angular average of the pair distribu- 
tion function of the hard-sphere sohd, /9^^^(ri, r2) [ssIIsgI] . 

We have followed the WCA prescription (after Weeks, 
Chandler, and Andersen Is?]) to split the potential 
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namely: 
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The properties of the reference system have been con- 
nected to those of the hard-sphere sohd through the "bhp 
function" formahsm [s^, i-e. the effective hard-sphere di- 
ameter (TwcA is determined by the imphcit relation: 



?/hs(r){exp[-/3i;pert(f')] -exp[-/3i;iis(r)]} dr = , (4) 



where Vhsir) is the hard-sphere potential (which depends 
on ctwca) and 2/hs(r) = 5hs(?') exp[/3Dhs(r)] is the corre- 
sponding cavity function. Another definition of the hard- 
core diameter (BH, ^^), namely: 



CTbh = / {1 - cxp[-/3i;pcrt(?')]} dr 

"'0 
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has been tested in this work. The expression 0, which 
corresponds to a first-order approximation to Eq. Q 34] , 
is appropriate to describe the rapidly rising repulsive in- 
teraction Wpcrt(?')- Other prescriptions for the separation 
of the potential (Q) (see Refs. and 0]) have been 
analyzed in this study; their predictions for the thermo- 
dynamic properties of the model at issue are, however, 
less accurate on the whole than those obtained through 
the WCA approach and will not be discussed further. 

In the results' section we also investigate a second- 
order correction to the free e nerg y, early proposed by 
Barker and Henderson in Ref. [4l|, and recently applied 
to systems with short-range interactions by Foffi and 
coworkers js^, namely: 
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where Phs is the pressure of the solid of hard spheres. 
The correction in Eq. ^ was obtained by dividing the 
space into concentric spherical shells, and assuming that 
the volume of each shell has the compressibility proper- 
ties of a macroscopic portion of the space (see Ref. 
for full details on the procedure and the approximations 
involved). 

We have resorted to several well-established simula- 
tion results in order to obtain the properties of the hard- 
sphere crystal. Specifically, we have used the Hall ana- 
lytical equation of state for the pressure 42], which ac- 
curately fits the molecular dynamics data of Alder and 
coworkers , ,43] : 
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where 77 is the packing fraction and ?7cp = 7r\/2/6 
is the close packing fraction; 7 = 4(1 — v/Vcp) and 
ao = 2.557 696, oi = 0.125 307 7, 02 = 0.176 239 3, 
03 = -1.053 308, 04 = 2.818 621, 05 = -2.921934, 
flg — 1.118 413. The free energy is obtained by ther- 
modynamic integration 
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where for the reference free energy, Fhs(p)j we have used 
the Frenkel and Ladd 'iij Monte Carlo result for the FCC 
crystal of hard spheres, namely: PF{^^/N = 6.537 9(09) 
at rj/rjc-p = 0.7778 in the thermodymamic limit. Their de- 
termination followed the Einstein crystal method, based 
on the construction of a reversible path from the solid 
under consideration to an Einstein crystal with the same 
crystallographic structure. As for the radial distribution 
function, 5hs(^)! we have used the analytical form, pro- 
posed by Weis [s^ and Kincaid and Weis 45], which 
expresses the distribution as a sum of gaussian peaks 
centered around the nearest-neighbour FCC lattice sites, 
with parameters calculated by Choi and coworkers (4^ 
to fit their simulation results. In Ref. i40| an analytical 
form for the cavity function yhsif) inside the hard core 
entering Eq. is also proposed, by directly extending 
to the solid the Henderson and Grundke scheme for the 
fluid phase [13 ]. 



C. Liquid state theories 

We have determined the thermodynamic and struc- 
tural properties of the fluid region through the 
MHNC m and SCOZA theories. We recaU that 
in the MHNC the Ornstein-Zernike equation, 



h(r) = c(r) + p / c{\v - T'\)h(r')dT' , 



(9) 



is coupled with the cluster expansion for the radial dis- 
tribution function g{r), namely: 



g{r) = exp[— /9?j(r) + h{r) — c{r) + B{r)] 



(10) 



where h{r) = [g{r) — 1] and c(r) are the pair and direct 
correlation functions respectively, and B[r) is the bridge 
function [23] ■ A closure to Eqs. Q and pO(l is obtained 
through the Percus-Yevick hard-sphere bridge function 
i3^^(r, cr^s), with the hard-core diameter crjj, fixed so to 
enforce the thermodynamic consistency of the theory (see 
e.g. [4^'). In particular we require the equality between 
the fluctuation pressure Pc and the virial pressure Pv, 
namely, 



— dp = Pv 

XT 



(11) 



In Eq. pi|l xo and xt are the ideal gas and the isother- 
mal compressibilities respectively; the latter is obtained 



4 



as the (7 = hmit of the Fourier transform of the di- 
rect correlation function, Xq/xt = [1 — pc{q = 0)], ac- 
cording to the fluctuation theory; the density integral is 
performed at a constant temperature. 

The consistency between and the pressure esti- 
mated via the energy route, Pc, can be alternatively en- 
forced along isochoric paths, namely: 
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In Eq. (|12|l the free energy is obtained by integrating the 
internal energy U with respect to the (inverse) tempera- 
ture at constant density, according to the relation: 
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where (p, T) is a thermodynamic state whose absolute 
free energy (3F{p, T) /N is known. In order to avoid to 
cross the liquid-vapor coexistence region, combinations 
of isothermal and isochoric paths are required to reach 
thermodynamic points located in the high-density, sub- 
critical region of the phase diagram. 

As shown in Ref . |2J| , the global consistency is in gen- 
eral more accurate than the local one; the latter proce- 
dure enforces the equality between the isothermal com- 
pressibility calculated via the fluctuation theory on the 
one hand, and via the virial route (i.e. by differentiat- 
ing Pv with respect to the density) on the other hand. 
The improvement of the global over the local consistency 
approach becomes crucial when the range of the particle 
interaction is very short, and/or at high densities and 
low temperatures; the free energy and the chemical po- 
tential in particular turn out to be almost quantitatively 
predicted 24]. 

As far as the SCOZA theory is concerned, this repre- 
sents an advanced liquid state theory that is based on 
a mean-spherical (MSA) type closure relation, replacing 
the prefactor /3 in the closure for the direct correlation 
function by a yet undetermined, state dependent function 
K{p, T); this function is fixed by the consistency require- 
ment between the ener gy a nd the compressibility route 
to thermodynamics [2^l27j . Thus for a hard-core poten- 
tial of diameter dhs and with an attractive tail w(r), the 
SCOZA closure relations to the OZ equation @ read: 
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In relation (|15l) c\is{r) represents the direct correla- 
tion function for the hard-core s yste m; we have used 
the Waisman parameterization [48j, i.e., Chs('") — 
i^o/r exp[— zo(r — (Xhs)] (with well-described density- 
dependent coefficients Kq and zq) which is known to 
give accurate results for the hard sphere system. The 
consistency requirement between the two thermodynamic 



routes mentioned above leads to the partial differential 
equation (PDE): 
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where Xrod = Xt/xo is the reduced (with respect to the 
ideal gas) dimensionless isothermal compressibility given 
by the compressiblity route and u is the excess (over ideal 
gas) internal energy per volume provided by the energy 
route. Solution of the above PDE leads to K{p, T) and 
thus to the structure and thermodynamics of the system. 

The SCOZA has proven to give reliable results for the 
location of the coexistence curve and stays — in con- 
trast to most liquid state theories — reliable even in the 
critical region (see for an overview). Many exam- 
ples for continuum systems as well as for spi n sy stems 
(in few cases including even binary mixtures 4!?|) have 
demonstrated its validity. Limiting our considerations on 
continuum systems, the SCOZA suffers from two draw- 
backs (at least in its present version), despite its many 
merits: as the formalism of the SCOZA largely benefits 
from the availability of the semianalytic solution of the 
MSA for a hard-core Yukawa system (including an ar- 
bitrary number of Yukawa tails) ,50], its application is 
limited to systems where the attractive part of the inter- 
atomic potential, w(r), can be approximated by a suit- 
able number of Yukawa tails, while the repulsive core has 
to be replaced by a hard-core potential. In the present 
work this has been done at the node of the potential, tq, 
as the SCOZA PDE becomes unstable for repulsive in- 
teractions. While in general the first problem does not 
represent a serious restriction, up to now no remedy has 
been found to include softness of the repulsive part of the 
interaction. 



III. RESULTS AND DISCUSSION 

In the first part of this section (Figs. I1I3|) we com- 
pare the perturbation theory for the solid phase with 
the corresponding Monte Carlo results Then, the 

PT free energies are combined with the simulation data 
for the fluid phase in order to determine the solid-fluid 
boundaries. In the second part fFigs. Bl^ we present a 
complementary approach where the MHNC and SCOZA 
theories are used to predict the liquid- vapor equilibrium, 
and — in combination with the simulation inputs for the 
solid phase — the full phase diagram. Hinging on 
both comparisons, we have calculated in the last part 
(Fig. O the phase diagram on the basis of purely theo- 
retical approaches. 

As far as the PT is concerned, the expression ||SJ) gives 
a hard-core equivalent diameter smoothly varying from 
CTbh = 0.970 nm at T = 2200 K to CTbh = 0.972 nm at 
T = 1800 K. Such a weak variation reflects clearly the 
stiff, rapidly varying nature of the repulsive part of the 
Girifalco potential. The more refined expression Q for 
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the hard-core diameter, (Twca, which adds a further de- 
pendence on the density, involves just a minor correction 
to (Tbh evaluated through Eq. 

We report in Figs. lll2l the PT free energy, chemical po- 
tential and equation of state along several isotherms in 
the solid phase, as obtained by using either CTwca or ctbh as 
hard-sphere diameter, together with approximations (j^J 
or © for the free energy. In the figures, the theoretical 
predictions are compared to the corresponding simulation 
results • As can be expected, the use of either ctwca or 
Ctbh produces similar results. The PT predictions gener- 
ally agree with simulation data all over the solid phase, 
especially in the chemical potential vs the pressure. The 
equation of state compares less favorably to Monte Carlo 
results (see Fig. 12), although wider discrepancies occur 
for p > 1.27 nm~^, i.e. almost outside the coexistence 
region, since the melting line is located between p ~ 1.25 
and p ~ 1.27 nm""^. Remarkably, the accuracy of the 
numerical estimates in the expression (0) results from 
the balance between two almost comparable — and rel- 
atively large with respect to their difference — contribu- 
tions of opposite signs. The PT predictions slightly get 
worse if the second order correction to the free energy 
(however small, since it hardly exceeds 2% of the main 
contributions) is taken into account. As observed in the 
original paper j39l |. this can depend on the semimacro- 
scopic derivation of the last term in Eq. ijHJl , that is more 
appropriate for longer interaction ranges, in such a way 
that a reasonably large number of particles fits into the 
region of attractive potential. The derivative of the pres- 
sure involved in Eq. ^ can as well represent a source of 
numerical uncertanties. 

We may anticipate that the fine reproduction of the 
chemical potential, as obtained through the first-order 
PT, eventually leads to the most accurate predictions 
for the phase behavior of the model. For such reason, we 
concentrate on the expansion (PI and show, in Fig. |3| the 
corresponding phase diagram; predictions corresponding 
both to CTwcA and to CTbh are reported. The coexistence 
lines are determined by using the simulation results of 
Ref. for the fluid phase. It appears that the PT gives 
accurate predictions for the coexisting temperatures and 
densities, with marginal effects related to distinct pre- 
scriptions for the hard-core diameter. By converse, a sys- 
tematic overestimate characterizes the coexistence pres- 
sure above the triple point (see bottom panel of Fig.|3Jl, 
especially if CTwca is employed. We argue in this latter 
case that the definition of the cavity function inside the 
hard core , directly borrowed from the theory of sim- 
ple fiuids ^46j. is to some extent unappropriate for the 
Ceo model; moreover, numerical uncertainties can occur, 
due to the more cumbersome procedure needed to evalu- 
ate CTwca through Eq. Q , which involves in particular the 
calculation of the slope of ^ijs(r) at close contact. 

In summary, it turns out that a PT first-order expan- 
sion of the model's free energy around a properly defined 
reference solid gives reliable predictions for the coexis- 
tence properties, and, more generally, for the solid phase 



behavior of the system at issue. This finding is notewor- 
thy in that it implies that the hard-sphere behavior dom- 
inates, to a large extent, the structure of the solid phase; 
conversely, as discussed in , the properties of the fluid 
phase are markedly affected by energetic aspects related 
to the attractive part of the potential. Our results fur- 
ther support the use of the PT to characterize the solid 
phase of systems interacting through short-range forces, 
already documented in Refs. [3lll33. l33l | . 

As far as the theoretical description of the ffuid phase 
is concerned, we report in Fig.|3|thc SCOZA and MHNC 
predictions for the free energy, the chemical potential, 
the pressure, and the internal energy, along with the cor- 
responding Monte Carlo results of Ref. [ll|. It appears 
that the MHNC predicts to a high accuracy the chemical 
potential as a function of the pressure, a quantity di- 
rectly related to the coexistence properties of the model. 
A satisfactory agreement with the simulation data also 
emerges for both the free and the internal energies, while 
a marked discrepancy affects the estimate of the pres- 
sure, especially in the high-density regime. On the other 
hand, as is visible in Fig. ^ the thermodynamic prop- 
erties of the model are generally overestimated in the 
SCOZA framework. 

The tendency of the MHNC to overestimate the pres- 
sure under the global consistency procedure, already ob- 
served in Ref. 24], appears somehow unexpected, espe- 
cially considering the overall good performances of this 
theory. We can conjecture that it might depend on an 
overestimate of the effective hard-sphere diameter ct^^, en- 
tering the bridge function as the adjustable parameter 
to enforce the thermodynamic consistency in Eq. I|ll() 
or (|12|l . Indeed, the value of the virial pressure de- 
pends on a delicate balance between contributions of op- 
posite signs. The use of the Verlet-Weis bridge functions 
(which fit available simulation data), instead of the an- 
alytical Percus-Yevick ones, might as well improve the 
predictions for the pressure. As for the SCOZA, we argue 
that the discrepancies with simulation data must be re- 
lated to the approximate treatment of the repulsive part 
of the potential. In particular, the substitution of the 
soft-core interaction with a purely hard-sphere potential 
for r < ro causes an enhanced repulsion (which reflects in 
the overestimate of the free energy and of the pressure), 
especially at high temperature, where shorter distances 
can be sampled by the system. We have indeed verified 
by few simulation runs along the isotherm T = 2100 K 
(not reported here) that the agreement between SCOZA 
and Monte Carlo pressures considerably improves, pro- 
vided that both approaches are compared on the same 
model, i.e. hard-spheres plus a Girifalco tail. 

We report in Fig.|31the SCOZA and MHNC predictions 
for the liquid-vapor coexistence properties of the model. 
The Gibbs Ensemble Monte Carlo (GEMC) binodal fine 
obtained in Ref. 15] with a sample of 1500 particles is 
also shown for comparison. It comes out that both the- 
ories reproduce quite faithfully the binodal curve, espe- 
cially as far as the vapor branch is concerned. The co- 
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TABLE I: MHNC and SCOZA critical and triple point den- 
sities (in nm~^) and temperatures (in K). Simulations: the 
critical point cor resp onds to the GEMC estimate with 1500 
particles of Ref. [isll : the triple point has been calculated in 
Ref. [H. 
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existing liquid densities are slightly underestimated, the 
SCOZA reverting this trend just below the critical point. 
The binodal curve obtained through the MHNC under 
the global consistency constraint definitely improves on 
the local one, as displayed in Fig. |S1 thus confirming the 
preliminary evidence reported in Ref. As for the 

SCOZA, an internal compensation between the overesti- 
mate of the chemical potential on one hand, and of the 
pressure on the other hand, could be at the origin of the 
fairly good agreement with the GEMC binodal curve. 

The MHNC solution algorithm converges to a ther- 
modynamic consistent solution up to T = 1900 K (see 
Fig. In view of early studies on HNC-type theories 
(see e.g. H^l), the difficulties to explore the critical region 
may be intrinsic to such approaches, rather than an arti- 
fact of numerical procedures. Consequently, the available 
MHNC binodal points must be fitted to some expected 
behavior in order to determine the critical parameters; as 
a common practice, we have used the scaling law for the 
coexisting densities to calculate the critical temperature, 

^cr' 

Pliquid - Pvapor CC \T ~ Tcrf , (17) 

where the exponent takes on the non-classical effective 
value P — 0.32. Once the critical temperature is known, 
we have calculated the critical density, pcr, by fitting the 
MHNC binodal points with the law of rectilinear diam- 
eter: Pliquid - Pvapor = Per + A\T - T„\, with A tO bc 

determined from the fit. At variance with the MHNC, 
within the SCOZA the critical region can be approached 
in principle with arbitrary precision and therefore nei- 
ther a fitting procedure, nor a scaling behavior hypothe- 
sis must be invoked to calculate the critical point param- 
eters. It emerges from Fig. |S1 that the theoretical pre- 
dictions closely bracket the GEMC critical point. The 
MHNC, SCOZA, and GEMC critical parameters are col- 
lected in Table n 

In Fig. we report the fluid-solid phase boundaries, as 
obtained by the combination of the MHNC and SCOZA 
predictions for the fluid phase and the simulation results 
of Ref. for the solid phase; for a more comprehensive 
overview, we also show a sketch of the binodal curves just 
discussed, as well as the full simulated fluid-solid coexis- 
tence points. As it is apparent, the overall good quality of 



the MHNC thermodynamic predictions provides a faith- 
ful reproduction of the phase boundaries of the model. A 
small overestimate of the fluid-solid coexistence pressure 
only occurs above the triple point temperature. As for 
the SCOZA, the solid-vapor equilibrium and the melt- 
ing line are satisfactorily predicted, while the observed 
shift of the chemical potential (see Fig. gives rise to 
a marked discrepancy in the coexistence pressure above 
the triple point; on the other hand, the overestimate of 
the pressure of the fluid phase is reflected in the shift of 
the theoretical freezing line toward lower densities. 

We have determined the MHNC and SCOZA triple 
points, also reported in Tabled from the intersection be- 
tween the corresponding binodal and freezing lines. It 
comes out that the MHNC predicts the existence of a 
stable liquid phase with a small shift of the triple and 
critical densities to lower values, with respect to the sim- 
ulation results. The SCOZA also predicts a liquid pocket 
for the model, although restricted to a narrower temper- 
ature and density range. 

It emerges that the MHNC theory for the fluid phase, 
and the flrst-order FT for the solid phase represent, 
among these envisaged here, the most accurate theoret- 
ical tools to enquire into the full phase portrait of the 
Ceo model. The ensuing phase diagram is reported in 
Fig. Ul The MHNC binodal line is directly drawn from 
Fig. |31 As is visible, the agreement with simulation data 
is semiquantitative; the comparison of Figs. El El and[7| 
demonstrates that the errors in the location of phase 
boundaries are essentially brought about by the MHNC, 
while the use of the FT for the solid phase does not intro- 
duce any appreciable discrepancy with simulation results. 
For such reason we do not expect a significant improve- 
ment of the overall appearance of the phase diagram if 
the SCOZA predictions for the fluid phase substitute the 
MHNC ones. In the bottom panel of Fig. [3 the over- 
estimate of the fluid-solid coexistence pressure magnifies 
slightly the common trend already observed for the FT 
and MHNC approaches separately. 



IV. SUMMARY AND CONCLUDING 
REMARKS 

We have presented a theoretical investigation of the 
thermodynamic properties and the phase behavior of 
the Girifalco Cgo model, in the framework provided by 
three refined theoretical tools, as the Modifled Hyper- 
netted Chain (MHNC) with a global thermodynamic 
consistency constraint and the Self-Consistent Ornstein- 
Zernike Approximation (SCOZA) for the fluid phase, and 
a Perturbation Theory (FT) with various degrees of re- 
finement for the free energy of the solid phase. 

The theoretical predictions are assessed against our re- 
cent simulation study 0. It turns out that the phase 
portrait of the model is predicted fairly well, in particular 
as for the solid-vapor coexistence and the melting line. 
By converse, the inaccuracies in the thermodynamic pre- 
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dictions sensitively affect the determination of the fluid- 
sohd coexistence pressure above the triple point. The 
MHNC and SCOZA liquid- vapor binodals are sufficiently 
accurate; the critical point estimates closely bracket the 
Gibbs Ensemble Monte Carlo datum |13 . In comparison 
with the simulation and the MHNC freezing lines (which 
turn to be practically superimposed), the SCOZA tends 
to underestimate the density of the liquid branch of the 
fluid-solid coexistence, with a corresponding shift of the 
triple point to lower densities and higher temperatures. 

We have interpreted our results in view of the ba- 
sic assumptions underlying each theory. It has emerged 
that the SCOZA mainly suffers from the substitution of 
the Girifalco soft-core repulsion at short distance with a 
purely hard-sphere interaction, which leads to a corre- 
sponding overestimate of the pressure and of the free en- 
ergy in the fluid phase. We have also conjectured that the 
MHNC pressure might improve if the Verlet-Weis bridge 
functions were adopted instead of the Percus-Yevick ones, 
as currently adopted in this work. As for the PT, it turns 
out that both a second order expansion of the free en- 
ergy, and a refined definition of the equivalent hard-core 
diameter, lead to a slight worsening of the theoretical 
predictions on the whole, which could be related to more 
cumbersome numerical procedures involved. 



We have eventually succeeded in predicting to a sat- 
isfactory degree of accuracy the phase diagram of the 
Girifalco model, on the basis of purely theoretical ap- 
proaches. We have selected for this purpose a sophisti- 
cated MHNC treatment of the fluid phase, and the flrst- 
ordcr PT for the solid phase. We expect that our con- 
clusions about the most suitable theoretical schemes to 
investigate the phase diagram of the Girifalco model, can 
be reasonably extended to other systems characterized 
by similar interparticle interaction laws. In particular, 
we plan to analyze the phase diagram of other higher- 
order C„>6o fuUerenes, combining both simulations and 
theoretical schemes. We have already produced several 
Monte Carlo results for the Girifalco C84 fuUerene, to 
be used for such purpose. The related calculations are 
currently in progress. 
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FIG. 1: PT free energy (top) and chemical potential (bottom) along three isotherms in the solid phase: first-order PT (Eq. @) 
with uwcA (Eq. (|1J, dotted line) or ctbh (Eq. full line), as equivalent hard-core diameters; second order PT (Eq. @) with 
awcA (dot-dashed line) or asH (dashed line). Symbols represent the corresponding Monte Carlo results |1L] . 



10 




FIG. 2: PT equation of state in the solid phase (hnes), displayed in order of decreasing temperature from top to bottom, 
according to the legend of Fig. Q Symbols: corresponding Monte Carlo results [Tll |. For clarity sake the isotherms T = 2100, 
1950, and 1800 K are displayed. 
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FIG. 3: First-order FT fluid-solid coexistence densities (top) and pressures (bottom) of the Ceo model, with ctwca (triangles) or 
(Tbh (circles) as hard-core diameters. Squares: Monte Carlo results H^. For completeness, the GEMC liquid- vapor coexistence 
densities (top) and pressures (bottom) are also shown as full lines p^ . The cross and the plus indicate the simulation critical 
and triple points, respectively. In the bottom panel dashed lines are guides to the eye for the Monte Carlo results. 
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Density [nm '] Density [nm" ] 

FIG. 4: MHNC (full lines) and SCOZA (dashed lines) thermodynamic properties along three isotherms in the fluid phase. 
Symbols: Monte Carlo results .11] at T = 2100 K (circles), 2000 K (squares) and 1900 K (diamonds). All curves are displayed 
from top to bottom in order of decreasing temperature. 
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FIG. 5: MHNC (dotted line with open circles) and SCOZA (dashed line) liquid-vapor coexistence; the GEMC results (full line 
with dots, S'lso shown. The crosses are the critical points. Pluses: MHNC under a local consistency constraint [l^ . 

The calculated coexistence points are shown as open circles (MHNC) and dots (GEMC) while the lines represent corresponding 
theoretical interpolations. 
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FIG. 6: MHNC (circles) and SCOZA (triangles) fluid-solid coexistence densities (top) and pressures (bottom). Squares: Monte 
Carlo results vTm with the corresponding triple point (plus). For completeness, the binodal curves drawn in F ig. El are also 
reported in the top panel. In the bottom panel the full line is the GEMC liquid-vapor coexistence pressure with the 
corresponding critical point (cross); dashed lines are guides to the eye for the simulation results. 
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theoretical fluid-solid coexistence densities (top) and pressures (bottom) of the Ceo rnodel. Squares: 



FIG. 7: Circles: fuUy 

simulation results vTrn . The MHNC and GEMC binodal curves with corresponding critical points shown in Fig. |3 are also 
displayed. In the bottom panel the full line is the GEMC liquid-vapor coexistence pressure (ISj with the corresponding critical 
point (cross); dashed lines are guides to the eye for the simulation results. 



